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We investigate the effect of strain-induced gauge fields on statistical distribution of energy levels of 
triangular graphene nanoflakes with zigzag edges. In the absence of strain fields but in the presence 
of weak potential disorder such systems were found in Ref. [1] to display the spectral statistics 
of the Gaussian unitary ensemble (GUE) due to the effective time-reversal (symplectic) symmetry 
breaking. Here show that, in the absence of disorder, strain fields may solely lead to spectral 
fluctuations of GUE providing a nanoflake is deformed such that all its geometric symmetries are 
broken. In a particular case when a single mirror symmetry is preserved the spectral statistics 
follow the Gaussian orthogonal ensemble (GOE) rather then GUE. The corresponding transitions 
to quantum chaos are rationalized by means of additive random-matrix models and the analogy 
between strain-induced gauge fields and real magnetic fields is discussed. 
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I. INTRODUCTION 

Graphene, a two-dimensional form of carbon, provides 
an intriguing condensed-matter analogue of (2-1-1) di- 
mensional quantum electrodynamics |2J. Several unique 
features of graphene, not found in conventional metals, 
semimetals, and insulators, follow from the effective de- 
scription for low-energy excitations given by the Dirac 
equation for spin- 1/2 fermions with zero rest mass [SHl]. 
Unlike a massless neutrino described by a similar Weyl 
equation [iBl , each effective quasiparticle in graphene has 
the electric charge — e and is coupled to external elec- 
tromagnetic fields via scalar and vector potential terms 
[7], offering several ways to tune the quantum states in 
order to control graphene-based electronic devices. Ad- 
ditionally, the dynamics of carriers in graphene can be 
affected by mechanical deformations inducing effective 
gauge fields [31 |H1 [H]- As strains exceeding 10% can 
be applied to graphene in a reversible manner [lOj . this 
new way of controlling the electronic structure has at- 
tracted significant theoretical [TTHTB] and experimental 
attention. Remarkably, the Landau quantization 
signalling the presence of strain-induced pseudomagnetic 
fields greater than 300 T was recently demonstrated [22] . 
For the opposite limit of weak strain fields, the existence 
of a zero magnetic field analogue of the Aharonov-Bohm 
effect in graphene [23] is predicted theoretically [16]. 

Surprisingly, the infiuence of strain fields on quantum 
chaotic behavior of electrons confined in graphene quan- 
tum dots [24J has not been discussed so far. Quantum 
chaotic behavior appears generically for systems, whose 
classical dynamics are chaotic, and manifests itself via 
the fact that energy levels show statistical fluctuations 
following those of Gaussian ensembles of random matri- 
ces [l5]. In particular, if such a system posses the time- 
reversal symmetry (TRS), its spectral statistics follow 
the Gaussian orthogonal ensemble (GOE). A system with 
TRS and half-integer spin has the symplectic symmetry 
and, in turn, shows spectral fluctuations of the Gaussian 



symplectic ensemble (GSE). If TRS is broken, as in the 
presence of nontrivial gauge fields, and the system has no 
other antiunitary symmetry j26| , spectral statistics follow 
the Gaussian unitary ensemble (GUE). For a particular 
case of massless spin-1/2 particles, it was pointed out 
by Berry and Mondragon |27j , that the confinement may 
break TRS in a persistent manner (i.e., even in the ab- 
sence of gauge fields) , leading to the spectral fluctuations 
of GUE [28]. 

When applying the above symmetry classification to 
graphene nanosystems [IK^S] one needs, however, to take 
into account that Dirac fermions in graphene appear in 
the two valleys, K and K', coupled by TRS If 
the valley pseudospin is conserved, a special (symplectic) 
time-reversal symmetry (STRS) becomes relevant, play- 
ing a role of an effective TRS in a single valley . Both 
real magnetic and strain-induced gauge fields may break 
STRS leading to the spectral fluctuations of GUE. As 
demonstrated numerically in Ref. fD, such fluctuations 
also appear for particular closed nanosystems in graphene 
in the presence of random scalar potentials slowly vary- 
ing on the scale of atomic separation. Such nanosys- 
tems include equilateral triangles with zigzag or Klein 
edges, i.e., with terminal atoms belonging to one sublat- 
tice. Generic graphene nanoflakes with irregular edges 
show spectral fluctuations of GOE [2S1 [31] , as strong in- 
tervalley scattering restores TRS in the absence of gauge 
fields. In contrast, the boundary effects are suppressed 
in open graphene systems, for which signatures of the 
symplectic symmetry class were reported |32|. 

In this paper, we analyze numerically the spectral 
statistics of triangular graphene nanoflakes with zigzag 
edges bent in-plane according to the strain geometry pro- 
posed by Guinea at al. (see Fig. [I]) . In the absence 
of strain flelds and the disorder, such systems were found 
in Ref. [J to show the Poissonian distribution of energy 
levels, which gradually evolves towards GUE when in- 
creasing the disorder strength indicating the transition 
to quantum chaos. Although the experimental energy 
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FIG. 1: Systems studied numerically in the paper. Bot- 
tom left: Triangular graphene nanoflake with zigzag edges 
characterized by the height H. Remaining plots: The same 
system bent in-plain employing the strain geometry proposed 
by Guinea et al. [T5] (top left) in the variant breaking all 
geometric symmetries (ac-strain), and (bottom right) in the 
variant preserving the mirror symmetry (zz-strain) . The radii 
of arcs limiting the flake area are R ± H/\^ for ac-strain or 
R± H/2 for zz-strain. The ratio H/R = 2 in both cases. 



resolution seems not sufficient as yet to discuss spectral 
statistics, regular graphene nanoflakes, including trian- 
gular ones bounded entirely with zigzag edges, were re- 
cently fabricated on metallic substrates |33j. Also, strain- 
induced gauge fields were shown to affect quantum states 
of similar nanosystems [22 . Here, we demonstrate nu- 
merically, that experimentally realistic strains may lead 
to clear signatures of quantum chaos even in a nanosys- 
tem which is otherwise perfect, i.e., not subjected to 
substrate-induced disorder, atomic-scale defects, etc. We 
also discuss the role of a particular strain geometry ap- 
plied, which may lead to the true or false STRS breaking, 
in a similar way as the geometry of chaotic Schrodinger 
systems leads to the true or false TRS breaking in the 
presence of magnetic fields [55] . 

The paper in organized as follows. In Secjllj we briefly 
discuss how magnetic fields and geometric strains con- 
tribute to the effective Dirac Hamiltonian for low-energy 
spectrum of graphene. In Sec. |IIH we present our main re- 
sults concerning the statistical distribution of energy lev- 
els of triangular graphene nanoflakes with zigzag edges. 
A quantitative comparison of the effect of uniform mag- 
netic flcld and the effects of different strain-induced fields 
on such a distribution is given in Sec. |IV[ The conclusions 
are given in Section \V\ 
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FIG. 2: Strain-induced gauge flelds appearing on the edges of 
triangular nanoflake for in-plane deformations of Fig. [l]in the 
limit oi R^ H (schematic). Red [blue] arrows are for Aj, > 
[Ax < 0]. Left: ac-strain leaves only to the approximate Sx 
symmetry, which is broken by a small Ay term; see Eq. (jsjl. 
Right: zz-strain breaks both Tsi and Sx, but it leaves the 
exact invariance under the reflection A — > —A combined 
with Sx, corresponding to the antiunitary symmetry To^ 
(false symplectic time-reversal symmetry breaking). 



II. STRAIN-INDUCED GAUGE FIELDS AND 
DIRAC FERMIONS IN GRAPHENE 

We start from the tight-binding Hamiltonian, includ- 
ing the nearest-neighbor hopping-matrix elements be- 
tween TT orbitals on a honeycomb lattice [31 [31] 
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where 5dij /do is the relative change in bond length, with 
do = a/\/3 being the equilibrium bond length defined 
via the lattice spacing in graphene a — 0.246 nm, and 
A is the vector potential related to the real magnetic 
field by B = rotA and incorporated as a Peierls phase 
(with the flux quantum $0 = h/e). Remaining param- 
eters are the equilibrium hopping integral to ~ 3eV 
and the dimensionless electron-phonon coupling /? — 
-d\og{t)/d\og{d)U=d, -2-2,^. 



A. The effective Hamiltonian and its symmetries 

It can be shown, that for low-energy excitations, low 
magnetic fields, and deformations slowly varying on the 
scale of atomic separation a, Htb Q reduces to the ef- 
fective Dirac Hamiltonian [H [HI [H] 



"Hcff = vpTocr ■ (p+eA) - vpTzCr ■ A, 



(2) 



where vp = 3todo/{2h) ~ lO^m/s is the energy- 
independent Fermi velocity, cr = (ax, cry), (Ji and 
(i — 1,2,3) are the Pauli matrices acting on sublattice 
and valley degrees of freedom (respectively), ctq (tq) de- 



notes the unit matrix, and p 



-iH(dx,dy) is the in- 



plane momentum operator. TieS ([2| is written in the 
so-called valley-isotropic representation [7|; namely, it 
acts on spinors -0 = [ipAjipB, —''P'bt ~'4^'a\'^ ' where A/B 
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is the the sublattice index, the primed and unprimed en- 
tries correspond to K and K' valleys. The strain-induced 
gauge fields can be written as 
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where c is a dimensionless coefficient of the order of unity, 
Uij = ^{diUj + djUi) = 1,2) is the strain tensor 
for in-plane deformations [3S], and we have chosen the 
coordinate system (x, y) such that the x axis corresponds 
to a zigzag direction of a honeycomb lattice. 

In the absence of gauge fields {A~A~Q), the Hamil- 
tonian ^ is invariant upon the antiunitary operations 

k 
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where C denotes complex conjugation. We notice that 
Ty = 1, and thus Ty represents the true time reversal 
coupling the two valleys; whereas 7^ — 7^ = —1, leading 
to the Kramer's degeneracy of the two valleys and to 
the additional Kramer's degeneracy in each valley. It is 
easy to see that real magnetic fields A 7^ break all 
the symmetries associated with 7fc-s. To the contrary, 
introducing the strain fields A 7^ and keeping A = 0, 
one may only break the invariance under the symplectic 
time reversal To (STRS), leaving TRS and the invariance 
under the valley exchange Tx unaffected f37^. 

The above discussion is complete, providing we ignore 
noncoUinear local magnetization (which may appear on 
the edges of a graphene nanoflake [3H1 EH] ) , so one can 
assume that the boundary condition to the effective Dirac 
equation "HcffA = Etp preserves TRS [7]. Recent first- 
principle study by Potasz et al. [59] shows, that local 
magnetization decays relatively fast when small graphene 
flake is charged out of the neutrality point. Although 
the extrapolation of this result onto much larger flakes 
as considered here may require some further analysis, it 
seems natural to expect, that the magnetization does not 
affect statistical properties of energy levels significantly 
[10] ■ Such an assumption if further supported by the fact, 
that no disambiguous effects of edge magnetization were 
observed experimentally even for the lowest-lying energy 
levels of graphene nanoflakes so far. For these reasons, 
the effects of possible edge magnetization are neglected 
in the remaining parts of the paper. 



with 9^^{x) = (V^x/H) arcsin [H/{y/3R}]. R is the 
bending radius of the deformation applied to equilateral 
triangle of the height H; see Fig. [T] The labels 'ac' and 
'zz' indicate that the strain following from Eq. ([5|, here- 
inafter referred as ac-strain, predominantly affects the 
bonds oriented along an armchair direction, whereas the 
strain following from Eq. ([6]) [zz-strain) predominantly 
affects the bonds oriented along a zigzag direction of 
a honeycomb lattice. The maximal strain in the triangle 
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for ac-strain, 
for zz-strain. 



(7) 



We notice, that ac-strain breaks all geometric symme- 
tries of the triangular nanoflake for any H/R > 0, while 
zz-strain preserves the symmetry with respect to a mirror 
reflection Sx ■ {x,y) {~x,y). Such symmetries of the 
deformations map onto the symmetries of effective gauge 
flelds A in the Hamiltonian |2| . Keeping only the terms 
of the order of R^^ and R^ one can express the strain 
flelds following from Eqs. ([5| and ^ ioi R> H as 
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In Fig. (2! we plot and on the flake edges only, as 
the flelds inside the flake smoothly interpolates between 
extreme values which are reached at the system bound- 
aries. It is clear that ac-strain breaks the symmetry with 
respect to Sx (depicted with dashed vertical lines) due to 
the term (Ay)ac oc R~^ . For zz-strain, Sx symmetry is 
also broken by the term {Ax)zz oc R^^. Both strain dis- 
tributions break STRS, as they are not invariant under 
the transformation A — >■ —A. However, A^^ exhibits the 
symmetry under the antiunitary operation ToSx (with 
[ToSx)"^ = — 1). In the absence of intervalley scattering 
and in the limit of quantum chaos, such a symmetry shall 
lead to false STRS breaking and spectral fluctuations fol- 
lowing GOE rather then GUE [21 [57]. 



III. 



SPECTRAL STATISTICS 



B. Deformations considered in the paper 

We further limit our discussion to the deformations 
earlier considered in Ref. [15) . namely 
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This section presents the central results of the pa- 
per, concerning statistical distribution of energy levels 
for Dirac fermions conflned in graphene nanoflakes, in 
the presence of gauge fields introduced in Sec. [IT] For 
the numerical illustration, we took an equilateral triangle 
with zigzag edges containing Nq = 32758 carbon atoms, 
corresponding to H — 270 do and the physical sample 
area of As — (29 nm)^. Such a system was previously 
found to show negligibly weak intervalley scattering in 
the absence of magnetic field, as well as spectral fluctua- 
tions following GUE in the presence of disorder smoothly 
varying on the length scale of a (see Ref. [Ij). 



4 



0.1 




0.0 



(a) 



2 

<!>/% 




FIG. 3: Evolution of energy levels for a triangular nanoflake 
containing A^c = 32758 carbon atoms with varying magnetic 
field (a), ac-strain (b), and zz-strain (c). 



A. The effects of magnetic and strain fields 

In Fig. [3j we plot energy levels of our model system ob- 
tained by numerical diagonalization of Htb as func- 
tions of the magnetic flux $ — AsB [Fig. Wa)], and 
the parameter H/R characterizing ac-strain or zz-strain 
[Figs. [Sj^b) orjsjjc)]. We limit ourselves to extreme cases 
of uniform magnetic field applied in the absence of geo- 
metric deformations, and two distinct deformations given 
by Eqs. (|5| and Q in the absence of magnetic fields, as 
intermediate situations simply combine the features of 
these extreme cases. First, it is clear from Fig. [sj^a), 
that uniform magnetic field splits the valley degeneracy 
of electronic levels [41j and leads to level crossings char- 
acteristic for integrable systems. To the contrary, both 
ac-strain and zz-strain preserve the valley degeneracy and 
lead to avoided crossings [see Figs. [Sj^b) and [Sj^c)] char- 
acteristic for chaotic quantum systems [5S] . A secondary 
difference between the effects of ac-strain and zz-strain 
on the electronic structure of the system considered is 
related to the fact, that different deformations lead to 
different pseudomagnetic fields. Namely, Eq. (Is]) leads to 
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with = if ^ = ac or 77^ = 1 if ^ = zz. It is clear 
from Eq. ([9|, that the strain-induced Landau quantiza- 
tion may be observed in small systems j42j for zz-strain 
only (see also Ref. [IS]). As this issue is beyond the 
scope of the paper, we only notice that the systematic 
drift of lowest-lying electronic levels towards the zero- 
energy Landau level is totally absent in Fig. [sj^b) and 
visible in Fig. [sj^c). For this reason, the arrangement 
of ac-strain seems particularly interesting when studying 
strained graphene nanosystems. In principle, ac-strain 
may allow one to discuss physical effects of strain-induced 
gauge fields in the absence of pseudomagnetic fields, in 
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FIG. 4: Integrated level-spacing distributions C'^^(S') for the 
same system as in Fig. |3] in the presence of uniform magnetic 
field (a)-(c) or ac-strain (d)-(f). The values of total flux <1> 
or the strain parameter H/R are specifled for each panel. 
Insets show nearest- neighbor spacings distributions P'^-'(S). 
Numerical results are shown with black solid lines. Remaining 
lines are for Poisson distribution with the degeneracy g = 2 
of each level [blue dashed], with g — 1 [blue dotted], and for 
the Wigner surmise for GUE [blue dash-dotted]. 



analogy to the Aharonov-Bohm effect appearing for real 
gauge fields in the absence of magnetic fields, providing 
the system geometry is modified appropriately |23| . 

The above-mentioned basic signatures of quantum 
chaos, accompanying strained-induced gauge fields, are 
further supported with spectral statistics presented in 
Figs.HHTl 



B. Level-spacing distributions 

First, we discuss the level-spacing distributions 
P'-^'^iS) and their integrals C^^^S) = P^^\S')dS' (see 
Fig.^, with P'-'^^S) (A: = 1, 2, ... ) being the probability 
that the quantity {p{E)){En+k — En)/k is located in the 
interval {S, S + dS). En+k — En is the distance between 
fc-th neighbors in the level sequence Ei ^ E2 ^ . . . , 
and {p{E)) is the average density of levels in the en- 
ergy interval {E,E + dE), which can be approximated 
by {p{E)) ~ n~^As\E\/{hvF)'^ for \E\ < to. The nu- 
merical results, shown with black solid lines in Fig. |4j 
are obtained for about 800 energy levels < i?„ < 0.5 tg 
|43) . The theoretical curves (blue lines) are given by 



Px{S), if.g = l, 

ld{S) + lPx{S/2), if.g-2, 
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FIG. 5: Second-neighbor level-spacing distributions P^^'(S') 
for the same system in same physical situations as in 
Fig. [4] Red solid lines show the best-fitted approximat- 



ing distributions Ppoi-2xPoi(A; 5") (i4l [panels (b), (c)] or 
-Ppoi-GUE(Afit; S") (f6l [panels (e), (TyTwith Aflt specified for 
each plot. Remaining lines are same as in Fig. [4] 
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FIG. 6: Level-spacing distributions P^^\S) {^-(c) and 
P^'^\S) (d)-(f) for the same system as in Figs. |3l|5| in the 
presence of zz-strain (with the parameter H/R specified for 
each panel). Numerical results are shown with black solid 
lines. Remaining lines are for Poisson distribution with the 
degeneracy <? = 2 of each level [blue dashed], the Wigner 
surmise for GOE [blue dotted], or for the best-fitted approxi- 
mating distributions Ppoi_GOE(Aflt; S) ( |15[ ) [red solid lines on 
panels (e), (f)] with Aflt specified for each plot. 



exp(-S'), 
PxiS)^{ (7r/2)5exp(-^^74), 
{32/n^)S^exp{-ASyTT) 



for Poisson, 

for GOE, (11) 

for GUE, 



where we have used the Wigner surmise approximating 
P^^\S) for the relevant ensemble of random matrices 
[25] . The evolution of spectral statistics C'-"'^-' (5') and 
P^^-'(S') with the magnetic field, illustrated in Figs.Ula)- 
(c) , confirms that the system energy spectrum gradually 
transforms from Poissonian with the twofold valley de- 
generacy (manifesting itself for <i> ~ 0) towards Poisso- 
nian without such a degeneracy (approached for <& ^ $o)j 
showing no signatures of quantum chaos in the absence 
of strain fields. The evolution of same statistics with the 
varying ac-strain and fixed $ = [Figs. |4|^d)-(f)] unveils 
the spectral fluctuations characteristic for GUE, starting 
from H/R > 0.1. 

The corresponding results for the second-neighbor 
spacing distributions P'^' (S) are presented in Fig.^ The 
theoretical expectations depicted with blue lines follow 



from Eq. ( 11 ) via 



2j^'dS'Px{2S-S')Px{S'), 
Px{S), 
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if .9 = 1, 
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(12) 

[For instance, Ppl{g=i{S) = 4S'exp(-2S').] Addition- 
ally, we have utilized additive random-matrix models of 
the form [3 gl] 



H{X) 



Hq + XV 



(13) 



allowing one to generate (for < A < oo) the spacing 
distributions interpolating between that of two distinct 
ensembles of random matrices. 



Three particular choices of Hq and V lead to: 
• A distribution interpolating between Poisson with g — 2 (reproduced for A = 0) and Poisson with g = 1 
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(reproduced for A = oo) [see Figs. [5]ja)-(c) 
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A distribution interpolating between Poisson and GOE [see Fig. [6] 
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A distribution interpolating between Poisson and GUE [see Figs. [5]jd)-(f)] 
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The coefficients a(A), u{X), and c(A) are chosen such that {S)x = /o°° SPx{X; S)dS — 1 for any value of A f4F; Iq{x) 
in Eq. ( 15 1 is the modified Bessel function of the first kind. 



Formally, Eq. ( 14 ) represents the exact expression for 
second-neighbor spacing distribution in case Hq is di- 
agonal random matrix with twofold degeneracy of each 
eigenvalue and V is diagonal random matrix without such 
a degeneracy. Similarly, Eqs. ( 15 ) and ( 16 ) represent the 



exact level-spacing distributions for 2x2 random matri- 
ces of the form H{X) ( [T3| , with Hq chosen as diagonal 
random matrix and V being a member of GOE or GUE 
(respectively). It was found numerically, however, that 
such distributions approximate, with a suprising accu- 
racy, the actual level-spacings distributions of large ran- 
dom matrices [2] as well as distributions obtained for 
various dynamic systems undergoing transitions between 
symmetry classes |25 l l44l l47] . 



The distribution -Ppoi-2xPoi(A; 5) (14 1, with the best- 
fitted parameter A = Agt , is capable of reproducing the 
evolution of P^-^^S) when the valley degeneracy is being 
split by external magnetic field in the absence of strain 
fields [see red lines in Figs.jsf^b) and[5jc)]. Similarly, the 
distribution Ppoi-GUE(Afit; -S*) ( flH] ) is capable of repro- 
ducing P'-^^S) when the system undergoes transition to 
quantum chaos induced by ac-strain [see red lines in Figs, 
[sje) and[5|f)]. In the latter case, P^'^\S) exhibits tran- 
sition Poisson-GUE analogous to the transition demon- 
strated numerically in Ref. [T] for the similar system with 
smooth potential disorder. 

The evolution of spacing distributions P^^''(S') and 
P^^-* (S) with zz-strain is illustrated in Fig. [6j This time, 
transition to quantum chaos manifests itself by a system- 
atic crossover of the actual spectral statistics obtained 
numerically (black solid lines) between the theoretical 
predictions for Poisson and GOE (blue lines) with the 
twofold valley degeneracy preserved. Also, the approxi- 
mating distribution Ppoi_GOE(Afit; 5*) ( [T5| is capable of 
reproducing P'^'^' (5) during the transition Poisson-GOE 
[see red lines in Figs.|6][e) and|6|f)]. 
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FIG. 7: Spectral rigidity A3(L) for the same system as in 
Figs. |3||6] in the presence of ac-strain (a) or zz-strain (b). 
Datapoints show the results obtained numerically for differ- 
ent values of the strain parameter HjR (specified for each 
dataset). Lines are the theoretical expectations for the rele- 
vant ensembles of random matrices. 



C. Spectral rigidity 

The fluctuations of more distant spacings between en- 
ergy levels of quantum system can be described in a com- 
pact way by the spectral rigidity (48) 



1 / t^l"^ \ 
(L) = -(Min / dx \M(xq+x) -ax- hf ), (17) 

L\(a,b) J _i^/2 ' 



where x = {M{E)) and N{E) denotes the number of 
energy levels such that < £"„ ^ E. Theoretical expec- 



tations Af°'^(L), 



Af°-)(L), 



and A3*^^^-'(L) are given 



explicitly in Refs. [TJdS]. For the case of f;-fold degener- 
acy of each energy level one finds immediately from Eq. 
(17) that the corresponding theoretical expression needs 



to be rescaled via A 



(^^^)(L)=g2Af)(L/5). 



The numerical values of ^z{L) for our model system 
are shown in Fig. [7] together with theoretical expectations 
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FIG. 8: Level-spacing distributions P^^\S) for the same sys- 
tem as in Figs. |3}|7|in the presence of uniform magnetic field 
(with the flux $ varied between the panels) and persistent ac- 
strain [panels (a)-(c)] or zz-strain [panels (d)-(f)], with the 
strain parameter fixed at H/R — 0.2 in both cases. Numeri- 
cal results are shown with black solid lines. Remaining lines 
are for GUE with the degeneracy g = 2 of each level [blue 
dashed], or for two independent GUEs [blue dotted]. 



D. Systems with persistent strains in external 
magnetic fields 

For a sake of completeness, we consider now the sys- 
tems of Fig. [l] with the strain parameter fixed at H/R = 
0.2 (the cases of ac-strain and zz-strain are studied sep- 
arately) placed in uniform magnetic fields varying in the 
range ^ $/$o ^ 4 [see Figs. [s] and |9]. Although the dis- 
cussion is still limited to the deformations given by Eqs. 
([5]) and ([g]), the universal nature of spectral fluctuations 
in the limit of quantum chaos allows us to believe that 
the effects which we demonstrate numerically in this sub- 
section may be observable for graphene nanoflakes with 
persistent (for instance, substrate-induced) strains, such 
as studied experimentally in Refs. |^ [55] . 

The numerical results for level-spacing distributions 
P^^'^ (S) are shown in Fig. [s] (black solid lines). The theo- 
retical predictions for GOE and GUE (with the valley de- 
generacy) are given by Eqs. ( |10|11 ) and drawn with blue 
dashed lines. The predictions for evel sequences follow- 
ing from two statistically-independent GOEs or GUEs 
(blue dotted lines) are given by expressions derived by 
Robnik and Berry [25] 
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FIG. 9: Same as Fig.[7]but in the presence uniform magnetic 
field ($/$o is specified for each dataset) and persistent ac- 
strain (a) or zz-strain (b); H/R = 0.2 in both cases. 



for the relevant ensembles of random matrices. Simi- 
larly as for the statistics P'^^'^^(iS') discussed above, the 
evolution of A3(L) with ac-strain clearly exhibits tran- 
sition Poisson-GUE [see Fig. [Tf^a)] whereas the evolution 
of A3(L) with zz-strain exhibits transition Poisson-GOE 
[see Fig. [7]jb)], with the valley degeneracy (g = 2) pre- 
served in the absence of magnetic field. 



Ex{S) 



l-erf(^V7r/2), 
cxp (-4S'V7r) 



for GOE, 



(19) 



-5* + erf {2S/y^) , for GUE, 



being the probability that interval S contains no energy 
level of a single sequence following GOE or GUE. The er- 
ror function erf(2;) = (2/-y/7r) exp{—t^)dt. Hereinafter, 
we suppose identical densities of energy levels in both 
components of a combined sequence. 

The datasets presented in Figs. [8]^a)-(c) illustrates 
a systematic crossover of the actual spectral statistics 
P^^^ (S) between the theoretical predictions for GUE with 
the valley degeneracy g = 2 and two independent GUEs 
driven by external magnetic field in the presence of per- 
sistent ac-strain. Similarly, a field-driven crossover be- 
tween GOE (with g = 2) and two independent GOEs 
is clearly visible in Figs. [8]jd)~(f) in the situation with 
persistent zz-strain. Additionally, both crossovers be- 
tween the relevant ensembles of random matrices are 
visualized with the spectral rigidity A^lL); see respec- 
tively Figs. [9f^a) and |9f^b) for the cases of persistent ac- 
strain and persistent zz-strain. [Notice that theoretical 
expectations for two independent GOEs or GUEs are 
""^^(L) = 2A^^^(L/2).] Remarkably, either P(i)(S') 
or A3(L) approaches the theoretical predictions for two 
independent GOEs (or GUEs) when $ > $0, and the 
(approximate) valley degeneracy no longer affects the 
electronic spectrum of the system. 

It is worth to stress that, apart from lifting up the 
valley degeneracy, weak magnetic fields do not alter the 
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FIG. 10: Least-squares fitted parameters Aflt for transitions 
between random ensembles demonstrated numerically in Sec. 



Ill [see Eqs. (14l-(16l], as functions of the magnetic fiux $ 
[solid sy mbo ls], or effe ctive fluxes "l>cff,^ [open symbols] defined 
by Eqs. (201 and 1 22 1. Dashed line represent the best-fitted 
power-low relation given by Eq. (211. Inset shows the data 



for strained nanosystems directly as functions of the strain 
parameter H/R, compared with the data obtained for real 
magnetic fields, in the log-log scale. 



symmetry classes of chaotic Dirac systems we consider 
here. In the absence of intervalley scattering, turning 
on the magnetic field simply transforms a system into 
a pair of two independent chaotic systems (one at each 
valley) each of which is showing the same symmetry class 
as the original system at zero field: Namely, the uni- 
tary if strain fields break STRS (the case of ac-strain 
of a generic nature) or the orthogonal if a mirror sym- 
metry leads to false STRS breaking (the special case of 
zz-strain). Chaotic graphene systems with strong in- 
tervalley scattering, earlier considered in Refs. I5T]. 
show standard transition GOE-GUE, exhibiting another 
striking difference between Dirac billiards and generic 
graphene flakes (with irregular edges). Nevertheless, the 
former still can be modelled effectively within particu- 
lar graphene nanoflakes, with terminal atoms belonging 
predominantly to one sublattice. 



IV. RELATION BETWEEN STRAIN FIELDS 
AND REAL MAGNETIC FIELDS 

We supplement our numerical study of quantum chaos 
in strained graphene nanoflakes by comparing, in a quan- 
titative manner, the effects of geometric deformations 
given by Eqs. ([5| and ^ with the effects of real mag- 
netic fields on the spacing distribution P^^\S). 

As discussed in Sec. |III[ a significant pseudomagnetic 
field per se, Bg oc H/R ([9|), appears for zz-strain only. 
However, as the actual spacing distributions p(^)(S') can 



be rationalized by Px{^{it,S) given (respectively) by Eqs. 
(14 1, (15), or (16) for the cases of real magnetic field. 



zz-strain, or ac-strain, the numerical comparison of the 
best-fitted parameters Afit (provided in Fig. 10) allows 
one to define the effective fluxes 



off, J 



$0 '4 ^^''^^ 



(20) 



for either the cases of ac-strain and zz-strain. The co- 
efficients are adjusted such that A^t for strained sys- 
tems [open symbols in Fig. [lO] follow the approximating 
power-law relation [red dashed lines] found for the case 
of real magnetic field 



Afit ^0.41(1) X [$/$o] 



1.53(7) 



(21) 



with the numerical values of parameters obtained via 
least-squares fitting (the standard deviation of a last digit 
are specified by numbers in parenthesis). This leads to 



10.5(2) and 



8.8(1). 



(22) 



We notice, that the dimensionless parameter c^z given 
by Eq. (22) corresponds to Czz$o/-4s — 43 T, what is 



numerically very close to the value reported by the second 
paper of Ref. [15 for the limit of Landau quantization. 
Therefore, the relation between uniform magnetic and 
strain-induced pseudomagnetic fields, following from the 
statistical description of the evolution of energy levels in 
weak fields by means of additive random-matrix models, 
appears to be consistent with the corresponding relation 
arising from transport properties in strong fields. 



V. CONCLUSIONS 

We have discussed the selected spectral statistics of 
triangular graphene nanofiakes with zigzag edges in the 
presence of strain-induced gauge fields. Such systems 
may show the spectral fluctuations following GUE of 
random matrices, providing the link between chaotic 
graphene flakes PI [5^ [31] and Dirac billiards for massless 
spin- 1/2 fermions [27, 28j. 

In the absence of disorder, strain fields associated with 
moderate in-plain deformations drive a highly-symmetric 
nanosystem into chaotic regime. Moreover, our results 
show, that the system symmetry class is related to the 
particular arrangement the deformation: In a generic 
case of the deformation breaking all geometric symme- 
tries the unitary symmetry class is observed, as the strain 
field also breaks the effective (symplectic) time rever- 
sal symmetry (STRS) in a single valley. To the con- 
trary, if a single mirror symmetry is preserved, we have 
only the false STRS breaking leading to the orthogo- 
nal symmetry class. (Such a physical situation has no 
analogue in graphene nanoflakes with substrate-induced 
disorder [1] or irregular edges [29l El], where all geo- 
metric symmetries are naturally broken.) It is worth 
to stress, that although particular strain arrangements 
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studied here are different from those in existing exper- 
iments [T7Vl22j . the nature of the results aUows one to 
expect, that the system symmetry class remains unitary 
(orthogonal) for an arbitrary strain arrangement leading 
to the true (false) STRS breaking, providing the deforma- 
tions are sufficiently small that the effective Dirac theory 
for low-energy excitations can be applied. 

In both cases of the unitary and the orthogonal symme- 
try classes, spectral statistics obtained numerically follow 
those of the relevant ensembles of random matrices (GUE 
or GOE), with the twofold valley degeneracy of each en- 
ergy level. When the real magnetic and strain fields are 
applied simultaneously, the valley degeneracy no longer 
applies, and our system displays spectral fluctuations fol- 
lowing two statistically-independent GUEs (or GOEs), 
one per each valley. 

The evolution of spectral statistics with weak strain 
fields, exhibiting the transition to quantum chaos, is ra- 
tionalized using additive random-matrix models. The 
functional relations between the best-fitted model param- 
eters and the deformation strength allow one to compare 
the effects of strain fields with the effects of real magnetic 
fields in a quantitative manner. The resulting numerical 



relation between uniform magnetic and strain-induced 
pseudomagnetic fields stays in agreement with the simi- 
lar relation, earlier obtained for strong fields leading to 
the appearance of Landau quantization (see the second 
paper of Ref. [T5]). 

Although our work was primarily motivated by the fab- 
rication of regular graphene nanoflakes |33) . it may also 
be possible to observe the effects which we describe in 
artificial graphenes, such as the arrays of GaAs/AlGaAs 
quantum wells |49j . or in the recently discussed acoustic 
analog of graphene [SD] ■ 
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